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MINIMUM  MEAN  SQUARE  ERROR  PREDT CT ION  OK  AUTOREGRESS 1 VE 
MOVING  AVERAGE  TIME  SERIES 

Ny  II.  Joseph  Newton  and  Marcello  Pagano 

Inrt i  tut  ■ '  .  <f  :  I'.r  t  : li.  ,  Ti'.rir  AW  Ibl  ’ t  u  and 
Popt..  of  Biontatir.tics,  Harvard  Uni  Density ,  U.d.A. 

Ei.'V.w r-dr. :  Autoregressive  Mov i rtg  Average  Time  Series;  ToepUtz  Matrix; 
Modified  Gholesky  Decomposition 

LANGUAGE 

ISO  Fortran 

DESCRIPTION  AND  PURPOSE 

Let  Y(l),  Y(T)  be  a  sample  realization  of  a  mixed  autoregressive 

moving  average  time  series  !Y(t),  t  =  0,  +  1,  ...}  of  order  (p ,q) (denoted 
ARMA( p,r()  ).  Thus 

p  q 

5  .t(j)Y(t-j)  =  [  p(k)t(t-k)  ,  t  =  0,  +  1,  ... 

j=0  k=0 

for  constants  p,q  ,  u(0)  =  B (0 )  “  1  ,  i(l),  ...,  n(p),  B(l),  ...,  B(q)  , 

where  t:  ( • )  is  a  white  noise  series  of  zero  mean  uncorrelated  random  variables 

2  VP  j 

having  variance  o  .  The  zeros  of  the  complex  polynomials  g(z)  =  )  _  pa(j)zJ 

rq  k 

and  h(z)  =  /^_^P(k)z  are  assumed  to  be  outside  the  unit  circle. 

Subroutine  MXPD  calculates  exact,  memory  h,  horizon  t,  minimum  mean 
square  error  linear  predictors  Y(t+h|t)  and  (optionally)  prediction  variances 
j  of  Y(t+h)  given  Y(l),  ....  Y(t)  fot  h  =  h^,  ....  h.,  and  t  *  t^,  ...,  t^. 
Thus  Y(t+h|t)  is  that  linear  combination  of  Y(l),  ...,  Y(t)  closest  to 

i  2 

Y(t+h)  in  the  mean  square  sense  and  o  ^  =  E{Y(t+h)  -Y(t+h|t)}  is  the 
attained  minimum  mean  square  prediction  variance. 

If  q  =  0,  Y  is  a  pure  autoregressive  process  of  order  p  (AR(p)),  while 
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if  p  =  0  ,  Y  is  a  pure  moving  average  process  of  order  q  (MA(q)).  Subroutines 
ARPD  and  MAPD  calculate  Y(t+h|t)  and  (optionally)  a"  for  the  ARip)  and 
MA(q)  cases  respectively.  Separate  subroutines  arp  used  for  these  cases 
due  to  significant  simplifications  in  their  algorithms  from  the  general 
A  RMA ( p , q )  case. 


NUMERICAL  METHOD 

The  method  is  based  on  numerous  special  p,  q.crlles  of  the  modified 

Cholesky  decomposition  (MCD,  see  Wilkinson  (1967))  of  the  autocovar i ance 

matrix  of  an  ARMA(p,q)  process.  (See  Newton  and  Pagano  (1930)  for  details), 

Let  [A].j  denote  the  (i,j)th  element  of  ,o  mac  c .  k  A  and  define 

A^  h  T0EPL{a(0),  ...,  a(K-l) }  to  be  the  (K*K)  symmetric  Toeplitz  matrix 

having  [A^J^  =  a(|i-j|)  .  Let  A^,  =  LKDKL^  ,  K  >  i  be  the  MCD  of  the 

symmetric  positive  definite  nested  matrices  A^,  A,  L._e.  [\_  ^  !  . 

=  ,  for  i,j  <  K  ,  Lr  is  a  unit  lower-  t'.auc...  (K*K)  matrix,  and 

D  is  a  diagonal  matrix.  Then  the  sequences  o;  R.n.r<ces  I,  .  D 
K-  K '  K.  \ 

and  are  also  nested  and  we  write  (for  example)  ae  (i,j)th  element 


of  any  L^,  for  K  >  max  (i,j)  as  [1,].^  . 


Thus  [D) 


11 


[A] 


11  ’ 


"-'.i  •  "*i,j 


i-r 

).  f  R  i  j  ^  i  o  i  p  p  r  i  •  i  „  • 


€-i 


i-1 


[Oli.  -  [A]u 


l  lD|f  f!l  I;  ,, 


1  •  K 


We  define  the  following  quantities  for  K  1 

a)  r  -  TOE  PL  { R  (0) ,  -  R..(K-D!  -  I...  if  „ 

/. ,  K  /.  /.  /  ,  K  /  .  K  ,  K 


Rj,(v)  -  E|7.(t)Z(t+v)  ) 


».  +  I.  . 


whe  re 
and  Z()  is  an  AR(p) 


-J- 


process  having  coefficients  n(l),  ...»  u(p)  and  white  noise 
variance  . 


b) 

II 

X  V 

(x(i),  .... 

x(k))t  =  i-;K 

Y  ,  where 

K 

T 

Y  = 
-K 

(Y(l),  .... 

Y  (K) )  . 

c) 

r 

X,K 

-  Elwl>  - 

T 

LX,KDX,K!iX,K  ■ 

Note  that  I'  „  = 

A  t  K 

d) 

II 

a*/ 

(e  ( 1)  ,  .... 

««))T .  l-;k 

*K  • 

-1 


,  -T 


Then 


Y ( t+h | t )  =  X(t+h|t)  -  l  m  (  j  )  Y (t+h-j | t) 

j  =  l 

2  **-1  2 

°t,h  k£0[LZLX]t+h,t+h-k^DX^t+h-k,t+h-k 


where  Y(t+h-j|t)  =  Y(t+h-j)  if  j  ^  h,  and 
X ( t+h | t )  = 


, t+h-k £ J;t+h_ktLz] t+h,£lLX]£, t+h-kj  | 


q 

lLx]t+h,t+h-ke(t+h"k)  •  h  “  1 . . 

k=h 


,  h  >  q 


Thus 


2 

to  obtain  Y(t+h|t)  and  o^  ^  for  h  =  h^,  ....  h,,  and  t  =  t^,  t 2 


,  -1 


we  need  4,t2+h2  *  4!t2+h2  *  LX,t2+h2  -and  DX,t2+h2  '  Significant  reductions 
In  computing  time  and  storage  requirements  in  obtaining  the  elements  of 
these  matrices  are  afforded  by  noting: 


NOTE  1  Computing  L„ 


-1 


z,t2+h2 


The  jth  row  of  l.'  is  given  by 

/. ,  K 


>'r 


( 1  ’  ®K-l)  ’ 


(u^jCj-l) ,  ....  aj_j(l),  1,  Qk_j)  »  2  £  j  <  p 


(0 


»  n(p).  •••.  1.  p+1  £  j  <  K 


(1) 


(2) 


-h- 


where  ^(£1  =  a(^)  if  k  >  p  and 

_  ni+i(,)  "  Vi(1+,)f>](,+l'1) 

u  .  (1 )  58  -  -  ~  ■  *  •  *  J  •  - - 

l~Vi(J+1) 


;  i  -  l . J^p 


Thus  there  are  only  p(p+l)/2  distinct  nonzero, iiononc  elements  in  1. 
K  >  p  +  1  and  X(l)  =  Y(l)  while 


-1 

Z,K  ’ 


X(j)  = 


j-i 


V(j)  +  l  a  (£)Y(.i-£)  ,  j 
£=1  J_1 


Y ( j )  +  l  a(£)Y(j-£) 
£=  1 


=  2,  ....  p 


(3) 


j  >  P 


NOTE 


2  Computing  h 

•  L  L 


Let  Y(0)  =  1,  y(1),  y(2),  ...  be  the  coefficients  of  the  MA(«>) 
representation  of  Z(*)  ,  i.e. 

j'1 

Y ( j )  =  -  j  a ( j-£) Y (£)  ,  j  >  1  . 

£=max(0, j-p) 


Then 


lLz]p+j,p+j-£  =  yii) 


i'A+i  k  =  '  I  «(r)[L  ] 

Zp+i'k  r=l  Z  p+J-r,£ 


,  0  <  £  <  j  <  K-p 

,  I  =  1 . P-1,  J  >  1 


I...  =  (L*)'1  ^  [1. 

Z ,  p  Z ,  p 


JZ]  .  J.  tVir^r.j-k’ 

}-l«  r=1-k-*-l  ,J 


P-1 


Lim 


1  im 

K - *  oo 


1  j  •» 


]-l 

K 


J,  ,L^ 


2  =  R7(0)/o2 


(4) 

(3) 

(6) 

k<j<p 

(7) 

(8) 

(9) 


-  >— 


y  Y  (k)  =  R_(0)/o^ 
k=0  L 


<.A0)  =  a2/  I!  (l-«2(j)) 

\i  =  l  J 


By  (5)  and  (6)  the  distinct  elements  of  I,_  are  contained  in  its 

Z,t.+h. 

T  L  L  T 

first  p  columns,  the  last  of  which  is  (0  1,  y(l),  ...»  y ) ) 

Thus  MXPD  finds  the  elements  of  L_  by: 

z,t2+h2 

i)  finding  L  =  (L~ 1  )_1  ,  ii)  calculating  Y(l) . Y(M.-p) 

7. ,  p  Z ,  p  1 

(by  (4))  where  M^-p  is  the  smallest  integer  such  that 


y  Y  (j)  -  R7(O)/02  <  6  ,  (12) 

J-0 

iii)  use  (6)  to  obtain  remaining  elements  of  rows  p+1,  of 

first  p-1  columns  of  L  ,  .  Note  that  (10)  says  that  such  an  M. 

''*t2+h2  1 

exists  (which  hopefully  is  much  less  than  (8)  ant*  (9)  say 

that  all  further  elements  of  the  first  p  columns  of  L„  ,,  are  arbitrarily 

Z,t2+h2 

close  to  zero.  Thus  there  are  M^p  elements  to  calculate  and  store  for 

LZ,t2+h2  ‘ 


NOTE  3 


Computing  L 


X,t2+h2  ,  X,t2+h2 


To  compute  Lv  .  and  Dv  .  ,,  we  note  that  (defining  an(0)  =  1) 

X,t2+h2  X,t2+h2  U 

[  r  ]  =  f  {  a,  .(j-m)  y  «,  .(i-OR-tf-1),  i,J>l 

J  I  m=max(l,j-P)  ^  £=max(l,i-p) 


R^di-jl) 


,  i,j  >  p,  |i-j  I  i  q 


,  1  f  j  1  P,  i  >  p,  and 
i  -  j  >  q  or  if 
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R  (v)  =  a  l  B(k)6(k+v)  ,  v  =  0,  . . . ,  q  .  (14) 

X  k=0 

Thus  the  distinct  elements  of  i'^  t  are  R^(0),  ....  R^q)  an<3 
elements  in  the  first  p  +  q  rows  and  p  columns.  Further  only  Ry(0),  ...» 
Ry(p+q)  are  required  to  obtain  rx  t  +h  .  These  elements  are  obtained  via 
subroutine  MXCV  by  solving  first  for  v  =  0,  ....  max(p,q) 

f  aUJR^v-j)  =  ?  B (k) Ry  (v-k)  =  0  ,  v  >  q  (15) 

j=0  -  k=v 

where  R^v)  =  E{Y(t)c  (t+v) }  =  6^0  if  v  >  0  ,  where  f>  is  the  Kronecker  delta. 


.  min(v,p) 

RYe(-v)  =  B(v)o  -  l  a(j)RYf;(-v+j)  ,  v  =  1,  ....  q  .  (16) 

j=l 

Then  Rx(max(p,q)+1) ,  Ry(p+q)  are  obtained  by  (15)  for  v  =  max(p,q)+l, 

....  p+q  . 

To  obtain  L  _  and  D  we  note  that  the  pattern  of  zeros  in 

A  j  1 2 +n  2  ^  1 2  '  n  2 

L„  .  is  the  same  as  that  in  the  lower  triangle  of  r„  and  that  th« 

Xjtj+f^  X,t2+n2 

required  elements  of  T  can  be  calculated  as  needed  (by  (13))  without 

X,  p+q 

storing  them. 

Thus  L„  .  is  calculated  and  stored  in  one  matrix.  To  obtain  the  q 
X,p+q 

nonzero  nonone,  elements  of  the  rows  of  the  rest  of  Lv  ,  we  have: 

a  >  1 2  ' ' *^2 


“■  [Vk,k-i  ■  «<J> 


*  j  ~  1 »  •  •  •  *  ^ 


Urn  [llx]K>K 


Thus  rows  p  +  q  +  1,  ...  are  calculated  and  stored  in  a  matrix  having  q 


columns  until  the  elements  of  L  ,  D  have  converged  (at  row  M_  say)  ,  i.e. 

XX  4 

|  SIX]M  ,  -  #(1-J«)I  '  ».  II0XI„  «  -  "'I  <  '•  J  "  V1 . V1  (17) 

2 ,  2  2 
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iv)  Find  and  f  ,jz  m  +p  ^  i  j  ’  1  '  1  i  )  1  P  by  (7),  (6), 

and  (12)  (stored  in  the  Integer  HONK  and  (M^*p)  matrix  AI.7.)  . 

v)  Find  R  (0),  ....  Rv(q)  by  (14)  via  subroutine  MACV  (stored  in 

X  X 

constant  RXO  and  q-vector  RX) . 

vi)  Find  L„  ,  and  I)v  ,  (stored  in  the  (p+q)*(p+q)  matrix  ALX1 
X,p+q  X,p+q 

and  in  the  first  (p+q)  elements  of  the  f^-vector  DX) . 

vii)  Find  M 2  and  [I,*x]  ,  1  =  P+q+1.  •••»  M2’ 

j  =  i-q,  ....  i- l  (stored  in  the  integer  MTWO  and  the  (M^xq) 

matrix  ALX2  and  the  rest  of  the  (F^-vector  DX) . 

viii)  Find  X  ,  e  by  (3)  and  (18)  (stored  in  the  t_-vectoisX  and  E)  . 
"t2  'l2 

ix)  Find  Y(t+h|t)  ,  for  h  =  h^,  ...»  hj  >  t  =  t^»  ...»  t2  (stored 

in  the  (t^tj+DO^-hj+D-vectors  YPD  where  Y(t+h|t)  =  YPD((t-t;[) 

(h2-h1^-l)+(h-h1+l))). 

2 

x)  Find  (optionally)  u  h  which  is  stored  like  YPD  in  a  (t^t^+l) 
(h2-hi+l)-vect°r  PVAR  . 


NOTE  6  Taking  advantage  of  convergence 

To  take  advantage  of  the  convergence  in  I,^»  and  the  user 

specifies  an  absolute  convergence  criterion  6  and  integers  IR0WS1, 

IROWS2  as  the  row  DIMENSION  of  arrays  ALZ  and  ALX2 ,  DX  respectively. 

Thus  IR0WS1  and  IROWS2  must  be  chosen  to  exceed  what  can  reasonably  be 
expected  to  be  M^  and  respectively  or  else  a  nonconvergence  failure 
indicator  will  be  returned  in  IFADi.T.  Of  course  if  1KOWSI  or  1ROWS2  are 
given  the  value  1 2  t',en  convergence  need  not  be  reached  for  the  algorithm 
to  finis  properly.  Setting  IROWSl  or  IROWS2  smaller  than  tjH^  allows  the 

possibility  of  obtaining  predictors  for  long  time  series  with  a  minimum 
amount  of  required  storage. 


-9- 


Further,  if  the  ^  are  not  to  he  calculated,  the  matrix  ALZ  is  not 
needed  and  l ROWS  1  can  be  set  equal  to  1. 

NOTE  7  Algorithm  for  ARPD 

For  q  =  0  we  have  r_  =  r  ,  r  =  D  =?  L  =  I  D 

Z.K  Y.K  X, K  Z,K  LX,K  *  UX,K 

=  K  '  T*1US 


Y  ( t+h  |  t )  «*  -  £  a  ( j  )  Y  ( t+h-j  1 1 ) 

j-1 


2  h_1  2 

a.  L  =  l  [L  ) 


t,h 


(DJ 


k=0 


Z  t+h, t+h-k  1  Z  t+h-k, t+h-k 


2  h_1  ? 

=  o  l  Y  (k)  if  t  >  p 

k=0 


NOTE  8  Algorithm  for  MAPD 


TOEPI 


For  p  =  0  we  have  1^  -  IR  *  Lz>|(  -  D.,^  -  IR  and  - 

(■VO).  •••»  ^(q) .  . ...  0)  ,  (so  that  calculation  of  ALX1  is 


avoided),  Ex  K  =  YK  •  Thus 


Y  ( t  +h 1 1 )  = 


Jh  [Lx]t+h, t+h-k  e(t+h-k) 


t  ,h 


B  (k)e(t+h-k) 

k=h 

0 

h-l 


l  IL 


k=0 


X ' t+h , t+h-k [ DX] t+h-k , t+h-k 


9  l'-1  ~ 

o2  l  B2  (k) 

k=0 


,  t+h  <  M„ 


,  t+h  >  M_ 


,  h  >  q 


,  t  <  M„ 


>  t  >  M 
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STRUCTURE 

SUBROUTINE  MXPD  (NP,  NQ,  ALPHA,  BETA,  S1CSQ,  Y,  I  OPT,  NT],  NT2 ,  NHI ,  N1I2,  NY PD, 
NPVAR ,  NPPNH2 ,  IROWS1,  IROWS2,  TROWS  1,  DEL,  RYE,  RYEO,  RY ,  RYO ,  IP,  YWK ,  RZO, 
RX,  RXO,  ALZI,  ALZ,  ALX1,  ALX2 ,  DX,  MONE,  MTWO,  X,  E,  YPD,  PVAR,  1 FAULT) . 

Formal  parameters 


NP 

Integer 

input : 

order  of  AR  part  of  model 

NQ 

Integer 

input : 

order  of  MA  part  of  model 

ALPHA 

Real  Array  (NP) 

i nput : 

coefficients  of  AR  part  of  model 

BETA 

Real  Array  (NQ) 

i nput : 

coefficients  of  HA  part  of  model  • 

SIGSQ 

Real 

input : 

variance  of  white  noise  in  model 

Y 

Real  Array  (NT2) 

input : 

data  vector 

I  OPT 

Integer 

input : 
1 

0 

option  switch  equal  to: 

if  both  predictors  and  variances 

to  be  calculated 

if  only  predictors  desired 

NT1 

Integer 

i nput : 

tj  (first  memory) 

NT2 

Integer 

input : 

t^  (last  memory) 

NHI 

Integer 

i nput : 

h  (first  horizon) 

NH2 

Integer 

input : 

h  (last  horizon) 

NYPD 

Integer. 

input : 

(NT2-NT1+1) (NH2-NH1+1 ) 

NPVAR 

Integer 

input : 

same  as  NYPD  if  IOPT  =  1,  >1  if  IOPT=0 

NPPNH2 

Integer 

input : 

NP  +  NH2 

I  ROWS 1 

Integer 

input : 

row  dimension  of  ALZ  in  calling  pro¬ 
gram  (^NP+2  if  IOPT  =1,  >1  if  IOPT=0) 
see  note  6  above 

I  ROWS 2 

Integer 

input : 

row  dimension  of  ALX2 ,  DX  in  calling 
program  ('NP+NQ-*1)  See  note  6  above. 

' ROWS  3 

Integer 

input : 

row  dimension  of  RY,  IP,  ALZI,  ALXI 
in  calling  program  (>NP+NQ) 

DEL 

Real 

input : 

absolute  convergence  criterion  (see 
(12)  and  (17)). 

RYE 

Rea]  Array  (NQ) 

output : 

(-D . Rye  (~q) 

RYEO 

Real 

output : 

^(0) 

RY 

Real  Array  (NP+NQ) 

output : 

Ry(  D . Ry(p+q) 

RYO 

Real 

out  put : 

Ry(0) 

!  C 

Integer  Array  (I  ROWS!) 

work : 

YWK 

Real  Array  (NPPNII2) 

work : 

ItZO 

Rea  1 

out  put  : 

Rj,(0) 

1  X 

Real  Array  (NQ) 

out  put  : 

Rx(l) . R/q) 

RXO 

Real 

out  put : 

Ry(0) 

M.Z  I 

Real  Ar ray ( I ROWS 3 , I  ROWS 3) 

out  put  : 

L'.p 

ALZ 

Real  Array  (I ROWS],  NP) 

output : 

if  IOPT  =  1,  ALZ  contains  [Ly)„, 

)  =  1 . P.  1  =  1 . , 

if  I  OPT  =  0,  ALZ  Is  not  used. 
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AI.X  1 

AI.X2 

Real  Army 

Real  Array 

( 1  ROWS  1 , 1  ROWS  J) 

( 1 KOWS2 , NQ) 

out  put  : 

out  put : 

'  X ,  p  tq 

|l-x  I  j  ^  <  1  =  p+n+i,  •  •  •  *  m2  > 

iiX 

Real  Array 

(1  ROWS 2) 

out  put : 

i  ”  <-<! . 1-1. 

f  I>x  1  f  .  i  =  1,  ....  M2 

Ml  INK 

Integer 

out  pul  : 

M  (see  note  6  above) 

m  rwo 

1 nt egor 

out  put. : 

M  (see  note  6  above) 

X 

Real  Array 

(NT 2)  ■ 

output : 

X(l) . x(t2) 

K 

Real  Array 

(NT2) 

output : 

e(l) . e(t2) 

YPD 

Real  Array 

(NYPD) 

output : 

( (Y (t+h | t) ,  h  =  h1,  ....  h2)  , 

• 

t  ~  tl*  ‘ - • *  t22 

PVAR 

Real  Array 

(NPVAR) 

output : 

if  IOPT  =  1,  ((o*  ,  h  =  h  , 

.  .  .  ,  h 2 )  ■,  t  t  j  ,  ....  t2) 

if  IOPT  =  0,  PVAR  not  used. 


1  FAULT  Integer 

output:  failure  indicator 

}■  Indi ea t i o n a 

Value  of  I  FAULT 

Meaning 

0. 

no  failure 

1 

NP<1,  NQ-'l,  or  IOPT  not  0  or  1 

2 

NTl-'NP+NQ  or  NT1>NT2 

3 

NH1<1  or  NH1>NH2 

4 

NYPD< (NT2-NT1+1) (NH2-NH1+1)  or  NPVAR<1 
or  IOPT  =  1  and  NPVAR< (NT2-NT1+1) 
(N1I2-NII1+1) 

5 

IROWS1-NP+2  and  IOPT  *  1  or  IR0WS1<1 
or  NPPNII2- NP+N112 

6 

I  ROWS2  •" NP4NQ+1 

7 

I  ROWS l^NP+NO 

a 

SIGSQ  0 

9 

Singular  matrix  in  Subroutine  MXCV 

10 

An  u  (J )  1  (see  (2)  ) 

11 

IOPT  =  1,  1 ROWR1 'NT2+NH2-NP,  and 
convergence  not  reaehed. 

12 

Nonpositive  |  I)  ]  .  .  encountered 

A  11 

13 

IROWS2^ NT2+NH2-NP-NQ  and  convergence 
not  reached. 

SUBROUTINE  ARPD  (NP, 

ALPHA,  SIGSQ, 

Y,  IOPT, 

NT1,  NT2 ,  NHl ,  NH2 ,  NYPD,  NPPNH2 

YWK ,  GAM,  YPD ,  PVAR, 

I  FAULT) . 

Formal 

pavamatars 

NP 

Integer 

input 

order  of  AR  model 

ALPHA 

Real  Array 

(NP) 

i  nput 

coefficients  of  AR  model 

SIGSQ 

Real 

Input 

variance  of  white  noise 

Y 

Real  Array 

(NT2) 

input 

data  vector 
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10PT 

Integer 

input : 

option  switch  equal  to 

1 

if  both  predictors  and  variances 

to  be  calculated 

0 

If  only  predictors  desired 

NT1 

Integer 

Input : 

tj  (first  memory) 

NT2 

Integer 

input : 

(last  memory) 

NH1 

Integer 

input : 

hj,  (first  horizon) 

NII2 

Integer 

input : 

h^  (last  horizon) 

NYPD 

Integer 

input : 

(NT2-NT1+1) (NH2-NH1+1) 

NPPNH2 

Integer 

input : 

NP  +  NH2 

YWK 

Real  Array  (NPPN112) 

workspace : 

gam 

Real  Array  (NPPNH2) 

output : 

y (1) ,  Y (NP+NH2) 

YPD 

Real  Array  (NYPD) 

output : 

( ( Y ( t  +h | t),  h  =  hj,  ....  h2). 

PVAR 

Real  Array  (NH2) 

output : 

if  I OPT  =  1,  a  ,  h  *  1,  . .  . ,  h 

L  ,  fl 

if  IOPT  =  0,  PVAR  not  used. 

I FAULT 

Integer 

output : 

failure  indicator 

Failure  Indications 

Value  of  I FAULT 

Meaning 

0 

no  failure 

1 

NP  <  1 

2 

NT1  <  NP  or  NT1  >  NT 2 

3 

NU1  <  1  or  NH1  >  NH2 

4 

SIGSQ  -  0 

5 

IOPT  not  0  or  1  or  NYPD 

<  (NT2-NT1+1) 

(NH2-NH1+1)  or  NPPN112 

<  NP+NH2 

SUBROUTINE  MAPD  (NQ, 

BETA,  SIGSQ,  Y 

,  IOPT, 

NT1 ,  NT2 ,  NH1 ,  NH2 , 

NYPD,  NPVAR, 

I  ROWS , 

DEL,  RX,  RXO , 

DX,  ALX,  MTWO , 

E,  YPD, 

PVAR,  I FAULT) 

ormal  parameters 

NQ 

Integer 

input : 

order  of  MA  model 

BETA 

Real  Array 

(NQ) 

input : 

coefficients  of  MA 

model 

SIGSQ 

Real 

input : 

variance  of  white  noise 

Y 

Real  Array 

(NT2; 

input : 

data  vector 

(OPT 

Integer 

input : 

option  switch  equal 

to: 

1 

if  both  predictors 

and  variances 

to  be  calculated 

0 

II  only  predictors 

desired 

NT1 

Integer 

input : 

tj  (first  memory) 

NT  2 

Integer 

input : 

t2  (last  memory) 

Nil  1 

Integer 

Input : 

h  j  ( first  hot .  zon) 

NH2 

Integer 

input : 

h2  (last  horizon) 

Nlll 

NH2 
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NY  Pit 

Integer 

input  : 

ril  .'-Ni  Ml)  (NI12-NII1+1 ) 

NPVAR 

Integer 

i nput : 

Same  as  NYPI)  i£  IOPT  =  1, 

■  1  if  IOPT  -  0. 

I  Rt)Wi> 

I nteger 

input : 

row  dimension  of  AEX  and  I)X  in 
calling  program  (see  notes  6  and 

H  above) 

DEL 

Rea  1 

input  : 

absolute  convergence  criterion 

RX 

Real  Array 

(NQ) 

out  pul  : 

R .,  ( 1  ) . Rx(q) 

RXO 

Real 

out  put  : 

A  v  ( ( ) ) 

I)X 

Real  Array 

(IROWS) 

out  put  : 

i  =  1 . m2 

A  l,X 

Real  Arrav 

( L  ROWS , NQ) 

out  put  : 

I'x’ij’  1  =  l’  M2’ 

j  =  i-q,  • • .  ,  i-1 

MTWO 

I nteger 

out  put : 

M  ,  (see  notes  6  and  8  above) 

F, 

Real  Arrav 

(NT2) 

nut  put : 

e  ( l ) . e  ( 1 2 ) 

YPD 

Real  Array 

(NYPD) 

output : 

( ( Y  (t+h  1 1 ) ,  h  =  h1 . h2), 

t  -  t .j  ,  1 2 ) 

PVAR 

Real  Array 

(NPVAR) 

output : 

if  IOPT  =  1,  ((o*  h),  h  =  hL . 

b^),  t  =  t^,  t2),  if  IOPT  =  0 

PVAR  not  used 

I  FAULT 

Integer 

output . 

failure  indicator 

b'ailure  Indiaationc, 

Value  of  IFAULT  Meaning 


0 

1 

2 

3 

4 


5 

6 

7 

8 


no  failure 

NQ  1  or  IOPT  not  0  or  1 
NT1  -  NQ+1  or  NT1  ">  NT2 
NH1  -  1  or  NH1  >  NH2 
NYPI)  <  (NT2-NT1+1)  (NH2-NH1+1) 
or  NPVAR  <  1  or  IOPT  =  1 
and  NPVAR  <  (NT2-NT1+1) (NH2-NH1+1) 

I  ROWS  <  2 
SIGSQ  <  0 

nonpositivo  [[)].,  encountered 

\  IX 

IROWS  <  MT2+NH2  and  convergence  not  reached 


A'.ixi Liary  af.<jorithmn 

SUBROUTINE  MACV  (NO,  BETA,  SIGSQ,  RX,  RXO,  IFAULT) 

rival  parame tern 


NQ 

BETA 

SIGSQ 


I nt cger 

Real  Array  (NQ) 
Real 


Input:  order  of  MA  model 

input  :  coefficients  of  MA  model 

input:  variance  of  white  noise 
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RX  Real  Array  (NQ) 

RXO  Real 

IFAULT  Integer 


output :  Rx(l) . ^(NQ) 

output:  R^,(0) 

output:  failure  indicator  equal  to: 
0  if  no  failure 
1  if  NQ  <  1 


SUBROUTINE  MXCV  (NP,  NQ,  M,  IROWS,  ALPHA,  BETA,  SIGSQ,  RYE,  RYEO,  WKM, 
IP,  RY,  RYO,  IFAULT) 


Formal 

parameters 

NP 

Integer 

input:  order  of  AR  part  of  model 

NQ 

Integer 

input:  order  of  MA  part  of  model 

M 

Integer 

input:  highest  lag  to  calculate 

(M  >  max  (NP.NQ) ) 

IROWS 

Integer 

input:  row  dimension  of  WKM,  IP 

in  calling  program 

(IROWS  >  max  (NP,NQ) ) 

ALPHA 

Real  Array 

(NP) 

input:  coefficients  of  AR  part  of  model 

BETA 

Real  Array 

(NQ) 

input:  coefficients  of  MA  part  of  model 

SIGSQ 

Real 

input:  variance  of  white  noise 

RYE 

Real  AYray 

(NQ) 

output:  RYe(-1) . RYc(-NQ) 

RYEO 

Real 

output:  Ry E (0) 

WKM 

Real  Array 

(IROWS, IROWS)  workspace: 

IP 

Integer  Array  (IROWS) 

workspace : 

RY 

Real  Array 

(M) 

output:  Ry(l),  ....  Ry(M) 

RYO 

Real 

output:  Ry(0) 

IFAULT 

Integer 

output:  failure  indicator 

Failure 

Indications 

Value  of  IFAULT 


Meaning 


0 

1 

2 

3 

4 


no  failure 

NP  '  1  or  NQ  <  1 

M  '  max  (NP,NQ) 

IROWS  <  max  (NP.NQ) 
singular  matrix  encountered 


The  subroutines  DECOMP  and  SOLV  as  described  by  Moler  (1972)  are  called 
by  subroutine  MXCV. 
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RESTRICTIONS,  I  I MK,  stouaok 

Tf  the  zeros  of  g(z)  are  not  outside  the  unit  circle,  then  one  of 
the  n  ( j )  will  be  greater  than  or  equal  to  one  in  magnitude  thus  giving 
IFAULT  =  10  in  MXPD.  If  the  zeros  of  h(z)  are  not  outside  the  unit  circle 
then  an  element  of  will  become  nonpositive  thus  giving  IFAULT  =  12  in 
MXPD  or  I FAULT  =  7  in  MAPD. 

The  bulk  of  storage  and  computing  time  in  MXPD  is  devoted  to  the 
'p  matrix  and  the  *  q  matrix  1.^.  The  values  and  increase 
as  the  smallest  zeros  of  g(z)  and  h(z)  approach  the  unit  circle. 
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SUUROUT  INE  MXPD! NP . NO. ALPHA .BETA.SIGSU.Y.IDPT.NT1.NT2.NH1 ,NH2, 
1 NYPO.NPVAR.NPPNH2.  I ROM  SI  ,1  ROMS 2. I  ROM S  3  .  DEL  . R  YE  •  R  YE 0.  R  Y.  R  VO »  IP • 
I  YMK  .RZO.RX.RXO . ALZ I .ALZ. ALX I .ALX2 . OX  *  NONE . NT WO • X .£ • VPD. P V AR • 

I  IF  AULT  I 

C 

C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  YPO  AND  (OPTIONALLY) 

C  PREDICTION  VARIANCES  PVAR  FOR  HORIZONS  NHI.....NH2  EACH  FOR 
C  MEMORIES  NTI.....NT2. 

C 

DIMENSION  ALPHA! NP ) . BE  TAINQ) . V<  NT  2) .R YE  I NQ I  »RVl  IR0MS3 ) . 

1  IP(  I  ROM S3  I  . ALZ I ( I  ROMS 3 , I  POM S3) , ALZC  (ROM  SI . NP ) , Y MK ( NPPNH2 )  • 

I ALX  I (  | ROMS  3. IR0MS3 ) . A L X2 ( I ROMS2 .N Q I  ,0X(  I ROM S2 ) . X ( NT2 ) *R X ( NQ)  , 

I  LI  NT 2  I . YPDI NYPD) ,P VARINPVAR ) 

DATA  ZERO. ONE#  EPS/ 0.0.  I . 0. 1 .E-l 0/ 

C 

C  CHECK  INPUT  PARAMETERS  : 

C 

I F  A  UL  T  =  l 

IF ( NP.LT. I . OR. NO.LT. 1 . OR.IOPT.LT. O.OR. IOPT.GT. I )  GO  TO  999 
If  A  UL  T  =  2 

IF ( NT  I .LE.NPFNQ.OR .NT2.L T.NTl )  GO  TO  999 
IF  AULT  =3 

IF<  NHt.LT.I. OR . NH2 . LT . NH 1 )  GO  TO  999 
IF  AULT  =A 

NCK  =  (NT2-NT  I  FT  )•( NM 2— NH | ♦ I  I 
IF( NVPD.LT.NCK)  GO  TO  999 
IF ( NPVAR.LT . 1 )  GO  TO  999 

IFINPVAR.LT.NCK.AND.IOPT.EQ.il  GO  TO  999 

I F  A  UL  T  =  5 

IF  I NPPNH2 .LI .NP»NH2 I  GO  TO  999 
IF<  IROMSI.LT. NPF2. AND. IOPT.EO. I  )  GO  TO  999 
IF  I  IROMSI .LT. I  )  GO  TO  999 
IFAULT  =6 

IE  I  IROMS2.LT. NPFNQ4-1)  GO  TO  999 
I  FAULT*? 

IF ( IROMS3.LT .NPFNQ )  GO  TO  999 
IFAULT  =  B 

!F< SIGSO.LE.ZERO)  GO  TO  999 

C 

C  FIND  RYO.RY! 1 ) .. .. .R Y(NPFNO)  Z 

C 

NPPNQ  =  NP  FNQ 

CALL  MXCVI NP.NQ.NPPNQ. I ROMS 3. ALPH A . BE T A . S I GSQ. H YE . R YEO. 

I ALX I . IP.RY.RYO.  IF!  ) 

I F  A  UL  T  =9 

IF(IFI.FO.A)  GO  TO  999 

r 

C  FIND  ALZI.RZO  (ALZI  INITIALIZED.  ALPHA(J.I)  FORMED  IN  ALXIIJ.I). 
C  J.I.LE.NP.  RZO  FORMED.  ALZI  FORMED  FROM  ELEMENTS  IN  ALX!)  Z 

C 

IF  AUL  T= l 0 
DU  20  I=I.NPPNQ 
DD  10  J- 1 . NPPNQ 
10  ALZ  I! I . J)  =  ZERO 
20  ALZ I ! I . I ) =ONE 
DO  30  I  -  1  .  NP 
30  ALX 1 ! NP. I ) x ALPHA! I ) 

C 

IF  I NP. F Q. I )  GO  TO  SO 
NP  M I -NP- 1 
DO  40  J=I.NPMI 


JJ=NPM| - J» l 
JJP 1  = JJ* 1 

PART=ALX1 ( JJPI  * JJP  I  ) 

IFf PART. GE. ONE)  GO  TO  999 
OEN  =ONE-PAR  T  *P AR  T 
OO  *0  I =1 . JJ 
JJPIMI=J JO|-I 

40  AL  XltJJ*tl=(  AL  XI  I  JJPI  •  I  1  — P  ART •  AL  X 1 ( J  JP I  »  J  J  P I  M I  )  1/ OEN 
C 

50  RZ  0=  S I GSQ 

OO  60  J=1.NP 

60  RZO=RZO/IONE-ALX1 t  J. J ) 4ALXI t J* J»  I 
C 

IF(NP.FO.l)  GO  TO  80 
DO  70  J=2.NP 
JM1=J-| 

OO  70  |=l,JM| 

J J=JMI-I Fl 

70  ALZItJ.llsALXMJMl.JJ) 

80  CONTINUE 
NPP 1 =NPF| 

OO  90  I=NPPl.NPPNO 
I FS  T  = I  —  NPP 1 

DO  90  J=1 .NP 
J J  =  IFSTAJ 
Jl=NP-JFl 

90  ALZ  It  I , J J 1  =  ALPHA ( J 1  I 

C 

C  IFflUPT.Ea.il.  FIND  MONE  ANO  ALZ  (INVERT  ALZ I .  FIND  ELEMENTS 
C  OF  NP  TM  COLUMN  OF  ALZ  UNTIL  CONVERGENCE. 

C  THEN  FILL  IN  REST  OF  ALZ I  : 

C 

IFt IOPT.EO.O )  GO  TO  230 
OO  110  1=1. NP 

00  100  J= 1 . NP 

lOO  ALZ( I . J)=ZERO 
110  ALZt  I  . I ) =ONE 

IF(NP.EO.l)  GO  TO  140 
OO  l JO  J=2.NP 
JMl=J-l 

OO  130  K  =  1 , JM 1 

C=ZERO 

JMK  = J-K 

JMK  P l  =  JMK ♦ 1 

DO  120  IR=JMKP1,J 

120  C=C-ALZt  J. IR 14ALZI < IR  *  JMK ) 

130  ALZ(J.JMKI=C 
140  CONTINUE 
C 

CK=R20/SIGSQ 
NPP  1=NPM 

AL  Z  t  NPP I , NP ) =  — ALPHA!  I  I 

SUM SO* ONE ♦ ALPHA ( I  1 • ALPHA  f 1  I 

MUP=  M I  NO  f IROVS1-NP.NT2FNH2-NP) 

00  160  J=2.MUP 
LLOtfsMAXOIO  .  J-NP1M 
LUP  =  J 
ALO  =  ZEHO 
CC=ONE 

00  150  LL=LLOW . LUP 


L=LL- I 


no  n  n 


i 


JML - J-L 
L 1 =  NP»L 

IFIL.GT.O)  CC=ALZ(L  1  .NPl 
150  AL 0= ALO— ALPHA ( J ML ) ACC 

NPP  JiKIP*  J 
AL7INPPJ.NP|=AL0 
SUMSQ=SUMSO»ALO*ALD 
IF  I ARSISUMSQ-CKl.LT.OCL)  GO  TO  13  0 
160  CONTINUE 

IF! MUP.LT.NT2FNH2-NPI  GO  TO  170 
MUNE =NT2FNH2 
GO  TO  190 
170  CONTINUE 
1FAULT=1 1 
MONE  = IROMS1 
GO  TO  999 
180  MONF=JFNP 
190  CONTINUE 
C 

IF(NP.EQ.l)  GO  TO  220 
NPM 1=NP— 1 
Ml MNP=M0NE-NP 
00  210  J  =  1  • Ml  MNP 
NPP  J— NPf-  J 

DO  210  L - 1 • NPM I 
C=ZCRO 

00  200  IR=l,NP 
NPPJMBsNPPJ- |R 

200  C=C~ALPHA( IR IAALZ1 NPPJMR.L) 

210  ALZINPPJ.L»=C 
220  CONTINUE 
C 

FIND  RXO.RXf  1  ) ... . »RX (NO  I  i 

CALL  MACVlNQ.BETAf  S I GS Q . RX • RXO • IF  1  ) 

FIND  ALX | .0X1  1  0X( NPFNOI  C  NOTE  THAT  ALPHAII.J)  FOR 

C  I  .LE.  J.LE.  I.LT.NP  IS  IN  ALZIIIM.I-jMI)  : 

C 

230  CONTINUE 

00  250  1=1 . NPPNQ 
OU  2  A  0  J  =  l, NPPNQ 
2 AO  ALX1I I  .J)=/ER0 
250  ALX 1  I  I . I l=ONE 
IF  A  UL  T  =  1  2 
C 

OX< 1 l=RY0 

0(1  3A0  t  =  2.  NPPNQ 

IMI =1-1 

OO  290  Jal.IMl 
C  *  ZERO 

00  260  L  =  l  .  I 
OO  260  M=1»J 
(LMM=I ABSIL-M) 

IF  (  ILMM.EQ.O  1  C<CMLZI(I  .L) AALZl I  J.  Ml  ARVO 
IF  I  I L  MM.  GT  .01  C  =  Cf  ALZIII.LIAALZK  J.MIARVl  ILMMl 
260  CONTINUE 

IF| J.EQ.l )  GO  TO  280 
JM I  *  J— I 

OO  2  70  L  = I • J  M 1 

C  =C—  ALX I ( I »L I  AO  X (L I A ALX 1 f  J • L ) 


270 


280  ALXI ( | . J ) =C/DX< J) 

290  CONTINUE 

C=RX0 

IF(l.GT.NP)  GO  TO  320 
C=ZERO 

OO  310  M=  1  ,  l 
C  I  =  ALZ I (  I  iM) 

C2-Z6H0 

OU  300  L  =  l.l 
ILMM=I ABStM-L) 

IF(  ILMM.EQ.O)  C  2  =  C  2F  ALZ I  (  I.LI4R70 
IFIILMM.GT.O)  C2=C2*ALZII  I,U*RY(  ILMM  I 
300  CUNTINUE 

310  C  =  C  FC I *C  2 

320  DO  330  L  =  t • I  Ml 

330  C=C-OX(L I  * ALXI C  t *L > * ALX 1  I  I . L ) 

IF(C.LT.EPS)  GO  TO  999 
340  OXIIMC 
C 

C  FINO  MTWO.  ALX2.  AND  THE  WEST  OF  DX  (IN  THIS  SECTION* 

C  I  AND  . . I-l  REPRESENT  THE  INDICES  OF  THE  NONZERO. 

C  NONONE  ELEMENTS  OF  THE  MATRIX  LSUBX.  NOTC  THAT  IF  INTEGERS 
C  M.LE.N.LE.NPFNQ*  THEN  LSUBX(N.M)  IS  STORED  IN  LSUBXIN.M) 

C  WHILE  IF  N.GT.NPFNQ.  THEN  LSUBXIN.M)  IS  IN  AL X2I N-NP-NO . 

C  NQ-(N-M)Fl),  M=N-NQ. ... .N- 1)  ' 

C 

IFAULT=12 

IUP=MINOt IROWS2-NP-NQ. NT 2FNH2-NP- NQ > 

DO  400  11=1. IMP 

I =NPPNQ F I  I 
I MNQ= I — NQ 

ALX2I  I  I  .  1  l=RXI  NO) /OX I  IMNO) 

IF(NQ.EQ.I)  GO  TO  370 
DO  360  JJ=2.NQ 
J= IMNQF J J-l 
NQl=NQ- JJFl 
J1 = J-NPPNO 
C=RXI NQI > 

JJMl =JJ- 1 

OO  350  LL-I.JJM1 
L= IMNQFLL-I 
J2=NQ-( J-LIF1 

IF1 J.LE.NPPNO)  C I = AL X  1 1 J , L I 
IFIJ.GT. NPPN  0 )  C1=ALX2< J1.J2) 

350  C=C-ALX2(  I  I.LLI40XID4CI 

360  AC  X2  I  I  I  • JJ)=C/DX(J  ) 

370  C=RXO 

OO  3B0  L=I.NQ 
LL* l-NQFL-l 

380  C=C-ALX2( II .L ) • ALX2 ( 1 I.L)«DXILL> 

OXII)=C 

IFIOXI I I.LT.EPSI  GO  10  999 

C 

DO  390  J J* I . NQ 
NQMJP1=N0-J JF| 

390  IFIAOSI  ALX2C  I  I  »  J  J  1 -BE  T  A  I  NOM  JPl  M.GT.DEL)  GO  TO  400 

IFI ABSIOXI  I  1 -SIGSO ) .GT.DEL )  GO  TO  400 
MTWO= I  I FNPPNQ 
GO  TO  420 
400  CONTINUE 

IF!  IUP.LT.NT2FNH2-NPPN0)  GO  TO  410 


n  n 


mtwo=nt2fnh2 

GO  TO  *20 
0  CONTINUE 
IF AULT= 13 
MT¥0= IHOWS2-NPPNQ 
GO  TO  999 
J  CONTINUE 

FIND  XII) . XI NT2) .E( 1 > .... .EINT2) 

NPP I =NP ♦ I 
XI  I  )  =  V«  1  I 

IFINP.EQ.il  GO  TO  *50 
on  *40  J=2 • NP 
C-YI J) 

JMl-J-l 

DO  *30  L=1.JM1 
JML  = J-L 

I  C  =  C  +  ALZ  I  I  J  *  JML  )  •  VI  J  ML  ) 

>  xi  j )=C 

>  OO  A/0  J=NPPl.NT2 
C=V( J) 

DO  *60  L=1.NP 
JML=J-L 

)  C  =C  *•  ALPHA  I  L  )  *  VI  JML  ) 

>  X(J|=C 

El  l  »  =  Xl  I  ) 

OO  *90  J=2 • NPPNQ 
c-=  x  (  j ) 

JM1 =J-1 

OO  480  L’li JM | 

>  C=C-ALXl< J.L)*E(L) 

>  FIJI  =C 
MI.OW-NPPNO*  1 
MUP=MINOINT2.MTWO) 

OO  510  JsMLOW.MUP 
JJ= J-NPPNQ 

C=  X  I J  ) 

OO  500  L=l .NO 
LL=J-NQ*L-I 

»  C  =C— ALX2 <JJ«L)*E  ICC  I 

>  E(J)=C 

IF( MUP.EQ.NT2)  GO  TO  540 
MUPP)=MUP»I 
OO  530  J*MUPPt .NT2 
C=X( Jl 

OO  520  L=l .NQ 
JML=J-L 

C=C-8£T  A ( L ) *E I JML I 

>  El J  )  =C 

i  CONTINUE 


ino  vpo  : 

NPOPT  sNH2— NH 1*1 
00  630  NT=NT|.NT2 
NSUF AH= I  NT -NT l  )*NPOPT 
NT  MNP  =  N  T  — NP 

OO  550  I  =  1  » NP 
I ( sNTMNP* I 


550 


T 

-  \ 

YWK (  I) =V<  I  I  » 

DO  610  NH=I.NH2 
NPPNH=NP*NH 
NTPNH=N I ♦ NH 
lROWLX=NTPNM-NPPNQ 
XTPH=ZERO 

IFINH.GT.NO)  GO  TO  590 
IFf NTPNH.GT.MTWO)  GO  TO  570 
00  56  0  K=NM.NQ 
I NDL  =  NQ—  K ♦ l 

inde=ntpnh-k 

X TPH=  XTPHFAL  X2  ( IROWLX.INDLI *E<  I NOE ) 

GO  TO  590 

DO  580  K=NH.NQ 
I NOE=NTPNH— K 

XTPH=XTPHF8E  T  A <  K ) *  E ( IN DE ) 

C*XTPH 

DO  600  J=l.NP 
I NOY-NPPNH- J 
C=C-ALPHA( J ) *YWK< I NOT  ) 

YWK(NPPNH) =C 
DO  620  l=NHl.NH2 
N 1 0=NP ♦  I 
N 1 l - I -NH l ♦ t 
N1 2=NS0F ARANI 1 
VPO ( N 1 2 ) =YWK(N10| 

630  CONTINUE 
C 

C  IF  IOPT.EQ.1.  RINO  PVAR  J 
C 

IF( IOPT.EO.O)  GO  TO  690 
Ml MNP= MONE  —  NP 
DO  680  NT -NT  I • NT2 
NSOF  AR= I NT-N  T 1  )*NPDPT 
OO  680  NM=NHI.NH2 
NINDX=NSQFARFNH-NH| M 
NTPNH=NT  FNH 
C=ZERO 

DO  670  KPl=I.NH 
K=KP1- 1 
KROW=NTPNH-K 
C 1=SIGSQ 

(FIKROW.LT. MTWO)  CI=0XCKR0W> 

C2=ONE 

IF (  K •  E Q.  0 )  GO  TO  670 
NPPK=NP»K 
C23ALZ(NPPK.NP» 

OO  660  I  R*  1 • K 
K  MR  =  K - 1 R 
C3=ZFPO 

IF (KMR.GT.M1 MNP I  GO  TO  640 
I N0L  =  NP  FK  MR 
C3=ALZ< INOC.NPI 
640  C4=ZEPO 

IFIIR.GT.NO)  GO  TO  660 
N10=Nf PNH-KMR 
IF  I N10.GT • MT  WO )  GO  TO  650 
NI I =NI 0-NPPNQ 
N1 2=NO-IH»I 
C  4  =  AL  X2  <  N 1 1.N12) 

GO  TO  660 


560 

5  70 

580 

590 

600 

610 

620 


6S0  C  4  =  BE  T  A  (  I  H  ) 

ft ft o  c?-C2*rj*c4 

670  t=C*CI»C2*C2 

680  PV  AH  <  N I  NOX  I  =  C 
C 

690  CONTINUE 

c 

IT  AUL  T=0 
999  HE TURN 
END 

SUBROUTINE  M  APD  <  NO  , BETA. SIGSO.Y. IOPT. NTI.NT2. NHI. NH2, 

I NY PD. NP VAR, IHOWS.DEL.RX.RXO rOX . AL X . MT MO ,E , YPD. P VAR • I F AULT ) 

C 


C  THIS  SUBROUTINE  CALCULATES  PREDICTORS  YPO  AND 
C  (OPTIONALLY!  PREDICTION  VARIANCES  PVAR  FOR  HOR12INS 
C  NHI . NH2  EACH  FOR  MEMORIES  NT1.....NT2. 


C 

DIMENSION  BETA(NQ) . Y(NT2  > .RX(NQ) . DXI  IRONS) • AL  X I  1 RO  MS*  NQ ) • 
!F( NT2 I . YPO< NY PD) .PVAR ( NP VAR ) 

OAT A  ZERO. ONE. EPS/0.0. 1*0.1 *E- l 0/ 

C 

IFAULT=I 

IF< NQ.LT . 1 .OR. IOPT .Lf* O.OR. 10PT.GT. I )  GO  TO  199 
IF  AULT =  2 

IF ( NT l .LT.NOFI .OR.NT1 .GT.NT2)  GO  TO  199 
IF AULT=3 

IF< NHI ,LT . I .OR. NHI .GT.NH2)  GO  TO  199 
IFAULT=4 

NC  K  =  (  NT  2—  N  T  1  ♦'I  )*(NH2-NHI*1  ) 

IF(NYPO.LT.NCK)  GO  TO  199 
IFINPVAR.LT. I)  GO  TO  199 

IF < NPVAR.LT.NCK. AND. IOPT .EQ.l )  GO  TO  199 
IF A  ULT  =5 

IF ( IRUWS.LT .2)  GO  TO  199 
IF  AULT=6 

IF( S I GSO.LE . ZERO )  GO  TO  199 
IF  AUL  T=  7 
C 

C  FIND  MTWO  ANO  ROWS  OF  AL X ,  DX  I 
C 

CALL  MACVINO.BETA. S I GS 0. RX • RXO • IF  I ) 

C 

OxC 1 )=RXO 
DO  10  1=1. NQ 
AL  X  (  I  .  I  >  =  QNE 
E<  I  1 =  Y (  I  ) 

IUP=MIN0( IR0WS.NT2 FNH2 ) 

DO  100  1=2, I UP 

| I sMAXOI I -NO- I .0)41 
IM 1=1—1 

NEL  TS= IMI-I l  F  l 

OO  30  J - 1 1 » 1 Ml 
J I  NO  = J— I  1+| 

1  M  J= I  - J 
JM | =  J— 1 
C=RX<  IMJ  ) 

IF  1 J.EO. It )  GO  TO  30 
Jl =MAX0I J-NO-I  ,  0 ) F  1 
DO  20  L= I  I  «  JM1 
LL=L-I 1  F  I 
J J=L  - J I ♦  I 


10 


n  n 


20  C =C-ALX (  I  ,LL  I *Dx <L )*ALX< J . J J > 

30  ALX1  l • J I  NO l—C /0 X ( J ) 

C  =  RXO 

on  40  J= 1  1  «  I  Ml 
J  I  N0=  J - I  I FI 

40  C=C-OX< J »*ALX( I ,JIND»**LX( (  .  Jf  ND> 

IF ( C  «LE «EPS I  GO  TO  199 
OX< I ) =C 


IF  ( 

I.GT. 

NT  2  ) 

GO 

TO 

60 

C=Y 

(  1  I 

DO  50 

J-  1 

.NELTS 

I  1  =  1- 

NFL  T 

SF  J 

-1 

50 

C=C-ALX<  I 

.  J) 

*61 

I  I  ) 

F  (  I 

)  =C 

I  F  ( 

I  «LE  • 

NQ) 

GO 

TO 

1  00 

60 

OO  70 

J=1 

•  NO 

• 

J  J=N0 

-  JF  1 

70 

IF  (  ARSC  AL 

X(  I 

.  J  ) 

-BETA! JJ  » ) • GE  *0  EL  1 

IF< A0S1OX1 I )-SIGSQ J.GE.DEL)  GO  TO  100 
Mr WO: I 

I F (  I • GE  *NT 2 )  GO  TO  110 
IP  1  =  I  ♦  1 

DO  90  J3 I P l • NT  2 

c=yc  ji 

DO  80  K*1.NQ 
JMK=J-K 

80  C=C-BETA|K)*EI JMK) 

90  EIJl-C 

GO  TO  120 
100  CONTINUE 
I F  A  UL  F  ■=  8 

IF( (UP.LT.NT2FNH2I  GO  TO  199 
110  MTWO-NT2>NH2 
120  1FAULT  =  0 

CALCULATE  PREDICTORS  J 

NPDPT-NM2-NH1 4  I 
OO  140  NT=NTl  t  NT2 
NSOF  ARs 1 NT -N  T 1 )*NPOPT 
DO  140  NH-NH1.NH2 
N I ND— NSOF  ARFNH—  NH1  F  t 
YPO|NINO)=ZERO 
IF(NH.GT.NQ)  GO  TO  140 
NTPNrt-NT FNM 
C=2ER0 

DO  130  K=NH,N0 
I NDE=NTPNH-K 
I  NDL  -  N  O—  K  ♦  l 
C  I  =8E  T  A I K  1 

IF! NTPNH.LE. MTWOI  Cl  =  ALX  < NT PNH* I NOLI 
C-CfC 1 4E (  INDE1 
YPD1NIND)=C 

IF  lUPTxl*  CALCULATE  VARIANCES  1 

IF( IOPr.EQ.0)  GO  TO  199 
DO  180  NT=NTI»NT2 
NSOF AR=(NT-NTI ) 4NPDPT 
DO  180  NM=NHI.NH2 
NHMI «NM- I 


130 

140 

C 

C 

c 


too 


n  n 


nind=nsofarfnh-nhi  *■  i 
N  T  PNH -N  f F NH 

c  =  s  ir.so 

IF  <  NT  »LT  •  MT  Mil  I  C^OMNTPNHl 
IF  I  NH.Ell.  I  )  GO  TO  t  80 
NHUP-MINOINHM1  .NO  I 
IF( NT. GE .  M  T  W  O  >  GO  TO  160 
DO  150  K-I.NHUP 
I NDD  =  NTP  NH-K 
INDL=NO-KU 

I  50  C  =  CFOX<  1NOO I FALX (NTPNH, 1NOL  > * AL X ( NTPNH • I NDL t 

GO  TO  ibo 

160  011  17  0  K=1.NHUP 

170  C-C F5 IGSQ*8E TAIK ) •  BE  T AIK ) 

lflO  PVAR(NlNOI=C 
C 

199  RFTUHN 

FND 

SUB ROUT INF  ARPDINP. ALPHA.SIGSQ.Y.  I  OP T , N T 1 . N T2 * NH 1 . NH2. 

I NYPD.NPPNH2 . YWK . GA  M«  YPO . P VAR. I  FAULT  > 

C 

C  THIS  SUBROUTING  CALCULATES  PREDICTORS  YPO  AND  < OP T 1 ON ALLY  I 
C  PREDICTION  VARIANCES  PVAR  FOR  HORIZONS  NH1.....NH2  EACH 
C  FOR  MEMORIES  NTI.....NT2 
C 

DIMENSION  Y  <  NT  2 1 • AL  PH  A ( NP I  . YKK ( NP  PNH2 1 • GAM 1 NPPNH2 )  . 

IYPOI NYPOl .PVAHINH? I 
OATA  ZERO.ONF/O.O. I .0/ 

C 

I F  A  UL  T  =  I 

IKNP.LT.tl  GO  TO  100 
IFAULT  =2 

IF ( NTI  .LE.NP. JR. NT  I .GT.NT2)  GO  TO  100 
IF AULT=J 

IF < NH1 .LT . I .OR.NHl .GT .NM2I  GO  TO  100 
IF  AULT-4 

1FI SIGSO.LE.ZERU)  GO  TO  100 
IF  AULT  =5 

IF(  10PT.LT. O.QR. IOPT.GT. 1 )  GO  TO  100 
I F  AUL  T  =  6 

NCK=<  NT2-NT I  FI  I *(NH2-NH1  *1  ) 

IF! NYPO.LT. NCK.OR.NPPNH2.LT. NPFNH2I  GO  TO  100 
I F  AUL  T  =  0 
C 

FIND  PREDICTIONS  : 

NPDPT  =NH2-NH I F I 
DO  50  NT  =NT l . NT  2 
NSOF AR=(NT-NT I IFNPOPT 
NT  MNP=NT-NP 

DO  10  1=1. NP 
I I=NT-NP ♦ 1 
10  YWK(  I  »  =  Y(  1  I  1 

DO  30  NH=I.NH2 
NPPNM=NP ♦ NH 
NTPNH=NT  FNH 
C -ZEPO 

DO  20  I-*l  ,NP 
( I =NPPNH- I 

20  C  =C  — ALPHA (l)AYWKCII) 

30  YWK  C  NPPNH 1 *C 


OO  *0  NH -  NH 1 • NH2 
I  NDNH=NSOFARFNH— NHI  M 
INOWK=NP>NH 

40  YPO! INDNHI -Y4K! iNDMKt 

50  CONTINUE 

IF  IOPT  =  l •  FINO  VARIANCES  : 

IF| IOPT.EQ.0)  GO  TO  100 
GAM (  1  ) =  —  AL  P  H  A  !  1  I 
IFlNH2.EQ.lt  GO  TO  80 
OO  TO  NH=2.NH? 

LLOM=MAXO!O.NN-NPt  F  | 

C-ZERO 

OO  60  LL=LLOW .  NH 
L  — LL— 1 
Cl -ONE 

IF(L.GT.O)  CI=GAM(LI 
NHML=NH-L 

60  C=C— ALPHA! NHML ) *C 1 

TO  GAM | NH ) =C 
80  PVAR! I l-SIGSQ 

IF  <  NH2.EQ. I  I  GO  TO  100 
OO  90  I  — 2  »  NH2 
IM  I  —  I  —  1 

90  PVAR! I  I— PVAR i I M I  I  F  S  I GS  Q*GAM ! I Ml J*GAM(  I Ml  I 
100  RETURN 
END 

SUBROUTINE  MAtvlNO.BETA.SIGSa.RV.RYO.IFAULTI 

this  SUBROUTINE  CALCULATES  MAINQ)  AUTOCOVARIANCES  OF  LAGS 
0.....NQ  INQ.GT.O) 

DIMENSION  BE TA! NQ I .RYl NQt 
DATA  ONE/ 1 . 0/ 

IFAULT=l 

IFINQ.LT.lt  GO  TO  40 
IFAUL  T  =  0 

C=ONE 

OO  10  l-I.NQ 
10  C=C fBE TAIII4BETA1I  ) 

RYO=C*SIGSQ 

DO  30  I  V  =  1  » NQ 
C  =  BET  A!  IV> 

IF! IV.EQ.NQI  GO  TO  30 
NQMI V=NU-IV 

OO  20  J*I.N0M|V 
JP I  V  - J ♦  I  V 

20  C-CFBETA! JI4BETA1 JPIV) 

30  RY< IV)=C*SIGSQ 
40  RETURN 
FNO 

SUBROUTINE  MXCV1NP .NQ. M, 1R 04 S . ALP HA . BE T A . S I GSQ . RYE . RVEO . 
I MKM . IP.RV.RYO. I FAULT! 


THIS  SUBROUTINE  CALCULATES  ARM A ! NP»  NO  t  AUTOCOVARIANCES  FOR 
LAGS  0 . • . • . M  I M.GE.MAX! NP.NQt »  NP.NQ.GT.Ot  i 


DIMENSION  ALPHA! NP ) . BET A 1NQ  )  .R YEI NO I .NKMI IRONS. IRONS)  • 
1 IP( IRONS) ,RV(M) 

OATA  ONE. ZERO/ 1 .0. 0.0/ 

C 

I  F  A  UL  T  =  t 

IF< NU.LT.i .OR.NP.LT. 1 )  GO  TO  110 
I F  A  UL  T  =  2 

MAXPQ=MAXO( NP. NQ) 

MM^MAXPQtl 

IF( M.LE.MAXPQ)  GO  TO  110 
IFAULT=3 

IF( IRONS.LT. MM)  GO  TO  110 
IF  A  UL  T  =* 

C 

C  FIND  HYE0.RYE1 l I . . • . . RYE (NO)  5 

C 

HYFO=Sl GSQ 
DO  30  tV*|.NQ 
C=S IGSQABETAIIV) 

NUP^M INOl IV  I NP ) 

DO  20  J= 1 .NUP 
I VMJ-I V- J 

IF( IVMJ.EQ.O )  GO  TO  10 
C=C-ALPHAI JiNRYEI IVMJ  ) 

GO  TO  20 

10  C  =  C—  ALPHA (J)*RYE0 

20  CONTINUE 

30  RYE(IV)=C 

C 

C  USE  DECOMP.  SOLV  TO  OBTAIN  RYO.RY I  1 )•••.. RYI MAXI NP . NQ ) ) 

C 

DU  40  I V= 1  «  MM 
RV (  IV  >=ZERO 

DO  AO  J - I . MM 
AO  NKMI IV. J)=ZERO 

C 

NPP |  =  NP  A 1 
NOPl=NO* 1 
DO  60  I VPl  =  l  *  NQP 1 
I  V  =  I  VP l- 1 
C=RYF0 

ir(IV.GT.O)  C  =CABE  T  A  f  IV) 

IF I  I V .EQ .NQ )  GO  TO  60 
DO  SO  K=IVPI«NQ 
KM I V  =  K  — I  V 

SO  C  =  C  ABE T  A  I K ) A RYE  I  KM I  V ) 

60  RY( I VPl )=C 
C 

DU  70  I  VP 1  =  1  .MM 
I  V* I VP l  —  l 

NKMI  I VPl . I VPl  ) =  NKM(  I VP  1.  1  VPl |  AONE 
DO  70  J=1.NP 
I  l-l ABSl  IV-JIA1 

70  NKMI  I  VPl  •  t  1  )MKM(  1  VPl.  I  1  )  A  ALPHA  I  J  ) 

C 

CALL  OECOMPI MM, IRONS. NKM.IP) 

1FI IP(MM).EQ.O)  GO  TO  110 
IFAULT  ~0 

CALL  SOLV I  MM  » I  RONS • NKM  «R Y* IP) 

RY  0=R V I  I  I 

on  80  I V*1 . MAXPQ 


IVP1=IVM 

80  RVI  IV)=Rr( IVPl  ) 

c 

C  USE  DIFFERENCE  EQUATION  TO  GET  THE  REST  OF  THE  RV  : 
C 

IFI M.6Q.MAXPQI  GO  TO  110 

DO  100  1 V=MM  »M 

C=2ER0 

00  90  J>1.NP 
I VM J=1V- J 

90  C=C— ALPHA ( J ) ♦R  Y |  IVMJ) 

100  RV (  IV  »=C 
C 


110  RETURN 
END 


END 

DATE 
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